*cd "C:\Users\eyurdagu\Dropbox\minW_insurance\codes\outputs\store\BM_June2022_FINAL"

use "pol.dta", clear
gen model="bm"

append using "..\BM_nominw_FINAL\pol.dta"
replace model="nominw" if model==""


append using "..\BM_nominwrecalx_FINAL\pol.dta"
replace model="nominw_recalx" if model==""

append using "..\BM_sgepsrecal_FINAL\pol.dta"
replace model="sgepscalib" if model==""


egen gpstate=group(a z tau ix)


gen tmodel=1 if model=="nominw"
replace tmodel=2 if model=="bm"

tsset gpstate tmodel

gen welfgain_nominw=100*((L1.val/val)^(1/(1-gamma))-1) if  model=="bm" //| model=="sgepscalib"

replace tmodel=0 if model=="nominw"
replace tmodel=1 if model=="nominw_recalx"


tsset gpstate tmodel

gen welfgainrecalx_nominw=100*((L1.val/val)^(1/(1-gamma))-1) if model=="bm" 


drop tmodel

egen gpmodel=group(model)

gen occ=1 if ix<=7
replace occ=2 if ix>7


*gen bite=0
*replace bite=1 if abs((log(wage)-log(minw))/log(minw))<0.0001


gen bite=1 if (wage-minw)/minw<0.1
replace bite=0 if (wage-minw)/minw>=0.1


gen meanval=.
forvalues  i=1/5 {
sum val [aw=distw] if gpmodel==`i'
replace meanval=r(mean) if gpmodel==`i'
}



gen mgain_nominw_ss=.
sum meanval if model=="nominw"
replace mgain_nominw_ss=100*((r(mean)/meanval)^(1/(1-gamma))-1) if model=="bm"


save "tempwithgains.dta", replace


gen welfgain_rel=.
gen welfgain_trunc=.
forvalues i=1/3 {
sum welfgain_nominw if gpmodel==`i' [aw=distw], d
replace welfgain_rel=welfgain_nominw-r(p50) if gpmodel==`i'

sum welfgain_nominw if gpmodel==`i' [aw=distw], d
replace welfgain_trunc=welfgain_nominw if gpmodel==`i' & welfgain_nominw>r(p1)  & welfgain_nominw<r(p99)
}

sum welfgain_nominw if model=="bm" [aw=distw], d
gen medwelfgain_nominw=r(p50) if model=="bm"


sum welfgain_trunc if model=="bm" [aw=distw], d

gen qa=.
gen qwage=.
gen qx=.

forvalues i=1/3 {
    
xtile qtemp = a if gpmodel==`i' [aw=distw] , nquantiles(10)
replace qa=qtemp if gpmodel==`i' 

xtile qtemp2 = wage if gpmodel==`i' [aw=distw] , nquantiles(10) 
replace qwage=qtemp2 if gpmodel==`i' 

xtile qtemp3 = x if gpmodel==`i'  [aw=distw], nquantiles(10) 
replace qx=qtemp3 if gpmodel==`i' 

drop qtemp qtemp2 qtemp3
	
}


gen ta=.
gen twage=.

forvalues i=1/3 {
    
xtile ttemp = a if gpmodel==`i'  [aw=distw], nquantiles(2)
replace ta=ttemp if gpmodel==`i' 

xtile ttemp2 = wage if gpmodel==`i' [aw=distw] , nquantiles(2) 
replace twage=ttemp2 if gpmodel==`i' 

drop ttemp ttemp2
	
}

egen gptemp=group(ix model) if model=="bm" //| model=="rho1"  //| model=="rho90"

gen taocc=.
sum gptemp
local ilast=r(max)
forvalues i=1/`ilast' {
xtile ttemp = a if gptemp==`i' & gptemp!=. [aw=distw], nquantiles(2)
replace taocc=ttemp if gptemp==`i' & gptemp!=. 
drop ttemp

}

gen daocc=.
sum gptemp
local ilast=r(max)
forvalues i=1/`ilast' {
xtile ttemp = a if gptemp==`i' & gptemp!=. [aw=distw], nquantiles(10)
replace daocc=ttemp if gptemp==`i' & gptemp!=. 
drop ttemp

}

/*
gen bitefirmbin=.
forvalues ii=0(0.2)0.8 {
replace bitefirmbin=`ii' if bitefirmbin==. & bitefirm<`ii'+0.0001
}

replace bitefirmbin=1 if bitefirm>0.80001*/

gen highbitefirm=(bitefirm>0.3)

gen welfgain_nominw_highbite=welfgain_nominw if highbitefirm==1
gen welfgainrecalx_nominw_highbite=welfgainrecalx_nominw if highbitefirm==1

gen welfgain_nominw_bc=welfgain_nominw if occ==1
gen welfgain_nominw_wc=welfgain_nominw if occ==2

gen welfgain_nominw_hibite_bc=welfgain_nominw if highbitefirm==1 & occ==1
gen welfgain_nominw_hibite_wc=welfgain_nominw if highbitefirm==1 & occ==2

preserve
collapse (median) occ welfgain_nominw welfgain_nominw_highbite logz logx [aw=distw] if (model=="bm"), by(ix daocc model)

gen ixocc=ix
replace ixocc=ix-7 if occ==2

forvalues cc=1/2 {


twoway (line welfgain_nominw ixocc if daocc==1, lcolor(red) lpattern(shortdash) sort) ///
(line welfgain_nominw ixocc if daocc==10, lcolor(blue) sort) ///
(connected welfgain_nominw_highbite ixocc if daocc==1, mcolor(red) lcolor(red) lpattern(shortdash) sort) ///
(connected welfgain_nominw_highbite ixocc if daocc==10, mcolor(blue) lcolor(blue) sort)  if  model=="bm"  & occ==`cc', ///
ytitle("Welfare gains (CE %)", size(large)) ylabel(, nogrid) xtitle("Skill group", size(large)) xlabel(1(1)7) ///
graphregion(fcolor(white) lcolor(white)) plotregion(margin(large))  ///
legend(order(1 "Low wealth, overall" 2  "High wealth, overall" 3 "Low wealth, high-bite firms" 4 "High wealth, high-bite firms") rows(2) region(lcolor(white))) 
graph export "comparefigs/welfgains_nominw_x_allvshighbite_occ`cc'_daocc.pdf", as(pdf) replace



}

restore



preserve
collapse (median) occ wage minw welfgain_nominw welfgain_rel welfgain_nominw_highbite  welfgainrecalx_nominw  welfgainrecalx_nominw_highbite ///
welfgain_nominw_bc welfgain_nominw_wc medwelfgain_nominw welfgain_nominw_hibite_bc welfgain_nominw_hibite_wc (p10) p10welfgain_nominw=welfgain_nominw (p90) p90welfgain_nominw=welfgain_nominw ///
[aw=distw] if model=="bm", by(ix model)

gen wagerel=log(wage)-log(minw)

gen ixocc=ix
replace ixocc=ix-7 if occ==2

gen welfgain_nominw_bc_norm=. if model=="bm" & ixocc==1
gen welfgain_nominw_wc_norm=. if model=="bm" & ixocc==1

gen welfgain_nominw_hibite_bc_norm=.
gen welfgain_nominw_hibite_wc_norm=.

sum welfgain_nominw_bc if model=="bm" & ixocc==1
replace welfgain_nominw_bc_norm=welfgain_nominw_bc-r(mean)

sum welfgain_nominw_wc if model=="bm" & ixocc==1
replace welfgain_nominw_wc_norm=welfgain_nominw_wc-r(mean)

sum welfgain_nominw_hibite_bc if model=="bm" & ixocc==1
replace welfgain_nominw_hibite_bc_norm=welfgain_nominw_hibite_bc-r(mean)

sum welfgain_nominw_hibite_wc if model=="bm" & ixocc==1
replace welfgain_nominw_hibite_wc_norm=welfgain_nominw_hibite_wc-r(mean)



forvalues cc=1/2 {
	
	
twoway (line welfgain_nominw_bc ixocc, lcolor(red) lpattern(solid) sort) ///
(line welfgain_nominw_wc ixocc, lcolor(green) lpattern(shortdash) sort)  if model=="bm", ///
ytitle("Welfare gains (CE %)", size(large)) ylabel(, nogrid) xtitle("Skill group", size(large)) xlabel(1(1)7) ///
graphregion(fcolor(white) lcolor(white)) plotregion(margin(large))  ///
legend(order(1 "Blue collar" 2 "White collar") rows(1) region(lcolor(white))) 
graph export "comparefigs/welfgains_nominw_x_allvsbcvswc.pdf", as(pdf) replace
	


twoway (line welfgain_nominw_highbite ixocc, lcolor(blue) sort) ///
(line welfgainrecalx_nominw_highbite ixocc, lcolor(red) lpattern(shortdash) sort) if model=="bm" & occ==`cc', ///
ytitle("Welfare gains (CE %)", size(large)) ylabel(, nogrid) xtitle("Skill group", size(large)) xlabel(1(1)7) ///
graphregion(fcolor(white) lcolor(white)) plotregion(margin(large))  ///
legend(order(1 "Counterfactual" 2 "Maintain wage levels") rows(1) region(lcolor(white))) 
graph export "comparefigs/welfgains_nominw_x_occ`cc'_allvsrecalx_highbite.pdf", as(pdf) replace



}


restore

